A single-domain spectral method for black hole puncture data 
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We calculate puncture initial data corresponding to both single and binary black hole solutions of 
the constraint equations by means of a pseudo-spectral method applied in a single spatial domain. 
Introducing appropriate coordinates, these methods exhibit rapid convergence of the conformal 
factor and lead to highly accurate solutions. As an application we investigate small mass ratios of 
binary black holes and compare these with the corresponding test mass limit that we obtain through 
a semi-analytical limiting procedure. In particular, we compare the binding energy of puncture data 
in this limit with that of a test particle in the Schwarzschild spacetime and find that it deviates 
by 50% from the Schwarzschild result at the innermost stable circular orbit of Schwarzschild, if the 
ADM mass at each puncture is used to define the local black hole masses. 
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I. INTRODUCTION 

The evolution problem of general relativity requires 
the specification of initial data that satisfies the Hamil- 
tonian and momentum constraints on the initial hyper- 
surface. There are different strategies to pose initial data 
for a specific physical situation, which typically involve 
a choice of free data and the subsequent numerical so- 
lution of the constraint equations to obtain the physical 
data; see pj for a recent review. An active area of re- 
search is concerned with initial data that describe two 
orbiting black holes 0, IM El IS El- F° r example, one 
can study the two-body problem of relativity by con- 
structing sequences of quasi-circular binary data sets that 
describe the quasi-adiabatic inspiral of two black holes 
HEUGIllIllEHH- Furthermore, binary black 
hole data sets are the starting point for evolutions in nu- 
merical relativity, e.g. [l5| . 

Important aspects of black hole data sets are the choice 
of hypersurface and how the physical singularity inside 
the black holes is treated. Concretely, since the con- 
straints give rise to elliptic equations, one has to specify 
a computational domain and boundary conditions. One 
possibility when considering two black holes is to work 
on R 3 with two balls excised. At the spherical excision 
boundary one can impose boundary conditions based on 
an isometry, as suggested by Misner [l^. This bound- 
ary condition is used in the first fully 3D numerical data 
sets and in the more recent thin sandwich type initial 
data sets ■ The excision boundary can also be defined 
by an apparent horizon boundary condition [TtI fl8| , see 
| id Il9| for recent applications. Other boundary condi- 
tions are motivated by Kerr-Schild coordinates [j]. 

Excising spheres introduces a technical complication 
into numerical methods on Cartesian grids. In finite dif- 
ferencing codes on Cartesian grids, the boundary points 
are not aligned with the grid and one has to construct 
appropriate stencils for a 'lego' sphere @, Hfl, l2i| • Alter- 
natively, one can work with adapted coordinates which 
match the spherical boundary, for example Cadez coor- 



or one can use multiple coordinate pat ches 
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dinates Q], 

wit h sp herical coordinates at the excision region 

An alternative to excision boundaries is to work on 
R 3 with two points (the 'punctures') excised, where the 
punctures represent the inner asymptotically flat infin- 
ity (Brill-Lindquist topology |23. l23j). Using the Brill- 
Lindquist topology directly is problematic numerically 
since one has to resolve a one over radius coordinate sin- 
gularity. However, it is possible to analytically compact- 
ify the inner asymptotically flat region, fill ing: in the miss- 
ing puncture points, and to work on R 3 [1, HI IS USUI- 
This simplifies the numerical method because no special 
inner boundary condition has to be considered 0. 

In this paper we focus on the construction of an effi- 
cient numerical method for the computation of black hole 
puncture data for vacuum spacetimes containing one or 
two black holes with linear momentum and spin. The 
numerical method, pseudo-spectral collocation, e.g. [28| . 
can give exponential convergence when the solution is in- 
finitely often differentiable (C°°). However, in its usual 
form puncture data is only C 2 at the punctures. We re- 
solve this issue by constructing an appropriate coordinate 
transformation that renders the puncture data smooth at 
the location of the punctures. Consequently, our pseudo- 
spectral method converges rapidly to highly accurate so- 
lutions, although the convergence rate is generally not 
exponential due to logarithmic terms in expansions at 
infinity, see below. 

Note that spectral methods have already been applied 
successfully to various elliptic problems in numerical rel- 
ativity, including neutron star initial data [29ll30ll3lll3^ . 
homo gen eous star models j33, 13; relativistic Dy son 
rings 351, and black hole initial data 
particular, 0,^1 use several coordinate patches to cover 
a binary black hole excision domain, and the situation 
can become quite complicated with 43 rectangular boxes 
and 3 spherical shells with various overlap and matching 
boundary conditions 0. While such multi-patch codes 
have a certain grid adaptivity built in (cmp. 'spectral el- 
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ements' |28|1. one of our motivations was to simplify the 
spectral method by construction of a simpler computa- 
tional domain. 

Therefore, a noteworthy feature of our spectral punc- 
ture method is that our choice of coordinates maps M 3 , 
including spatial infinity and the two puncture points, to 
a single rectangular coordinate patch. This is the case 
for spherical coordinates with a radial compactification, 
but recall that in addition we want to ensure smoothness 
at the punctures. 

The paper is organized as follows. In Sec. [n] we de- 
scribe our spectral method for the solution of the Hamil- 
tonian constraint on a single domain. After we introduce 
the puncture data in Sec. IIIII Sec. IIVI discusses analyti- 
cal issues and numerical results for a single puncture. In 
Sec. we develop our single-domain spectral method for 
two punctures and present the key result, which is rapid 
convergence of this scheme to highly accurate solutions. 
The application of our method to the case of small mass 
ratios and a comparison with a semi-analytic test mass 
limit can be found in Sec. IVII Finally, in Sec. IVIII we 
compute binding energies of puncture data in this limit. 
We conclude in Sec. IVIIII 
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Hence, the grid points Ai and Bj are the zeros 
the Chebyshev polynomials T nA (l — 2x) and T nB {— x), 
respectively, whereas the Lp k represent the zeros of 
sm(n v tp). Our spectral expansion is thus a Chebyshev 
expansion with respect to the coordinates A and B, and 
a Fourier expansion with respect to p. 

The spectral method enables us to calculate first and 
second derivatives of U from the values f/yfc at the above 
grid points within the chosen approximation order which 
is given by the numbers (n A , tib, n v ). Thus, for a vector 



II. THE SPECTRAL METHOD 

As will be discussed in detail in the subsequent sec- 
tions, for both the single and the two-puncture initial 
data problems an elliptic equation of the form 

f(u) = Au + g(u) = (1) 

arises for a function u. Here, A denotes the Laplace 
operator, and g is a source term which in general depends 
on u. 

In what follows we will introduce coordinates (A, B, ip) 
with 

Ae[0,l], Be [-1,1], ^G[0,27r), (2) 

for each specific case that we consider, in which u is well- 
defined within the spatial domain, in particular at its 
boundaries. 

As will become clear below, u always obeys a physical 
fall-off condition at spatial infinity, 

lim u = 0. (3) 

r — >oo 

Since in all cases to be considered, the coordinate A is 
introduced such that 

r — > 00 A — ► 1, (4) 

we consider an additional function U which is given by 

u=(A-l)U. (5) 

In our spectral method, the values of this function U 
are calculated at the grid points (Ai, Bj, ipk), i-e. 

U ijk = U(A i ,B j ,tp k ), (6) 



U — (U 000 , ■ ■ ■ ,U(n A -l)(n B -l)(n^-l)) T (H) 

we may fill another vector 

f(U) = (/oOOj • • • j f{n A -l){n B -l){n v -l)) T (12) 

by the evaluation of f(u) at the grid points (Ai, Bj, (fk)- 
This results in a non-linear set of simultaneous equations 

10) = (13) 

for the unknown Uijk- 

For the function U we find particular boundary con- 
straints by considering the elliptic equation ft), written 
in terms of (A, B, tp) at A = 0, A = 1, B = ±1. A solu- 
tion U that is regular with respect to A, B and (p must 
obey these requirements, which therefore replace bound- 
ary conditions that usually need to be imposed. These 
boundary constraints are called 'behavioral' [28| . In addi- 
tion, a desired 27r-periodicity with respect to <p is already 
'built-in' by the particular choice of our basis functions. 
Accordingly, using the spectral method with the above 
interior collocation points, no further work with respect 
to the boundaries needs to be done, for regularity and 
periodicity will be realized automatically. Hence, there 
is no other requirement constraining the function U. It 
is uniquely determined by the elliptic equation ft). 

For the numerical solution of the discrete equivalent, 
Eq. (|13|) . we address its non-linearity by performing 
Newton-Raphson iterations. The solution U is written 
as 

U = lim U N , (14) 

W^oo 

Un+i = Un-Vn, (15) 
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where Vn satisfies the linear problem 

JnVn = b N (16) 

with 

df - 

J N = -4,(U N ), b N = f(U N ). (17) 
oil 

There are different ways to solve the linear system aris- 
ing from multi-dimensional spectral methods, although 
some effort has to be made to obtain an efficient method 
since the one-dimensional spectral differentiation matri- 
ces are not sparse and the conditioning of the system can 
be problematic; sec for example [lj|,|28j. We solve (|lfc>|) 
with the preconditioned 'Biconjugate Gradient Stabilized 
(BICSTAB)' method (36|, and the choice of precondi- 
tioner is crucial for the overall efficiency of the method. 

We construct a preconditioner which is based on a sec- 
ond order finite difference representation of J/y. To this 
end we consider the linearized differential equation cor- 
responding to on the equidistant grid in coordinates 
(a, f3, ip) with 

,4 = sin 2 a, B = -cos/3. (18) 

Apart from the uniform distribution of our grid points, 
these coordinates have the additional advantage that U 
becomes symmetric with respect to the planes a = 0, 
a = 7r/2, p = and (3 = tt. Therefore, it is possible to 
calculate second order finite differencing approximations 
of first and second derivatives at any grid point by taking 
into account adjacent neighboring points only. 

The resulting matrix has at most seven non- vanishing 
entries per row and is therefore well suited for the applica- 
tion of a sparse system solver. We use the program pack- 
age ' hyp re' which offers a variety of sparse matrix meth- 
ods [33. A choice that works well in this context is the 
'Generalized Minimal Residual (GMRES)' method pre- 
conditioned with the algebraic multigrid code 'Boomer- 
AMG' [13. 

Our implementation of the above procedure uses the 
BAM code as infrastructure [l5j|. Although both BAM 
and hypre support parallelization, we have not paral- 
lelized the elliptic solves for our spectral method. Com- 
putation of the binary black hole solution shown in Fig.[S] 
takes four minutes on a Xeon/Linux workstation. 

In what follows we evaluate the convergence of the 
spectral method by computing the 'global relative ac- 
curacy' defined by 

S n ,m(U) = max(^ )B , v )|l - U n /U m \, (19) 

where U n denotes a specific nth order spectral approx- 
imation of the function U. The maximum is typically 
evaluated over a regular grid of 6 3 points. We take 
iia = Ub = 2n v = n, with the only exception n v = 4 
which we use in axisymmetric situations. 

Choosing a large value m defines a reference solution to 
which solutions at a lower order of approximation n < m 



can be compared. This method gives a reliable charac- 
terization of convergence in our examples. Furthermore, 
we reduce the error in the solution of the discrete non- 
linear system l|13[) below the error due to the finite order 
of the spectral approximation which therefore dominates 
the accuracy of the method. 

The convergence rate of a spectral method is called 
exponential if the logarithm of the total error of an ap- 
proximate solution depends linearly on the corresponding 
approximation order for sufficiently large order. This be- 
havior is usually encountered if the underlying solution to 
be approximated is analytic everywhere on the spectral 
domain. However, if the solution is only C fc -differentiable, 
the logarithm of the total error depends linearly on the 
logarithm of the approximation order. In particular, the 
slope of this line is [k + 2), and the scheme is called al- 
gebraically convergent to (k + 2)-th-order. Hence, from 
the numerical convergence of the spectral method one 
can deduce the differentiability of the solution to be ap- 
proximated (see Figs. 1, 2 and 4 below for representative 
examples corresponding to puncture initial data). 



III. PUNCTURE DATA 

In the ADM- formulation of a '3+1 '-splitting of the 
spacetime manifold, the vacuum Hamiltonian and the 
momentum constraint equations of general relativity read 
as follows: 

/. ,2 • 1\~ l\,J\ ,J = 0, (20) 
V(K ij - j ij K) = 0. (21) 

Here 7^ is the 3-metric, the extrinsic curvature, K 
its trace, and R, V are the Ricci scalar and the covariant 
derivative, respectively, associated with 7^. 

Following York's conformal-transverse-traceless de- 
composition method pj, we make the following assump- 
tions for the metric and the extrinsic curvature (8ij de- 
notes the three-dimensional Kronecker symbol): 

7ij = (22) 
Kij = 4>~ 2 + V itj - I 6 i:i divV^j . (23) 

The initial data described by this method are conformally 
flat and maximally sliced, K = 0. With this ansatz the 
Hamiltonian constraint yields an equation for the confor- 
mal factor ip, 

A?p + - ip 5 K i:j K ij = 0, (24) 
8 

while the momentum constraint yields an equation for 
the vector potential V, 

AV + - grad(divV) = 0. (25) 



4 



One can proceed by choosing a non-trivial analytic so- 
lution of the Bowen-York type for the momentum con- 
straint, 



IV. SINGLE-PUNCTURE INITIAL DATA 

For a single puncture at the origin of a Cartesian grid, 
we introduce spherical coordinates (r, ip) via 



N p 

V = J2 

n=l 



4av 



X n *P % 



x n x S r 



with poles at a finite number of N p spatial points, the 
locations of the punctures. Here the vector parameters 
P n and S n can be identified with the physical linear 
and angular momenta of the nth puncture. The vector 
x n points from the nth puncture to the point (x, y, z) 

•En {x X n7 1/ Vn^ 

norm. 



x = r cos 

y = r sin i? cos ip, 

z = rsin^sinyj, 



where 



rG[0,oo), i?g[0,tt], (pe[0,2n). 



z n ) T , and \x n \ is its Euclidian The conformal factor for a single puncture is 



In Q it is pointed out that for the extrinsic curvature 
determined by (|26[) a particular solution of the Hamil- 
tonian constraint is obtained by writing the conformal 
factor ip as a sum of a singular term and a finite correc- 
tion u, 



ip = 1 



N p 

\ - m r . 



u, 



(27) 



with u — > as \x n \ — > oo. The parameter m n is called 
the bare mass of the nth puncture. 

The main point of the puncture construction is that in 
terms of u the Hamiltonian constraint becomes a well- 
defined equation on the entire Cartesian 3-space (see [11] 
for a general existence theorem for such asymptotically 
flat initial data). However, it turns out that u is in gen- 
eral only C 2 at the punctures, although it is C°° elsewhere. 

As discussed in Section^ such a drop of differentiabil- 
ity implies that a spectral method can only be expected 
to be algebraically convergent to fourth order. We first 
show for a single puncture that a simple coordinate trans- 
formation can resolve the differentiability problem at the 
location of the punctures. After that we discuss similar 
techniques for two punctures. 

Note that by virtue of Theorem 1 by Dain and 
Friedrich [3^1 ■> the conformal factor can only be expected 
to be globally C°°-differentiable with respect to our co- 
ordinates (A,B,tp) if the individual linear momenta P n 
vanish. In fact, we will find that for punctures with lin- 
ear momenta the conformal factor possesses logarithmic 
terms when expanded at infinity, i.e. at A = 1. This holds 
also true if the total linear momentum, i.e. the sum of all 
P n vanishes. Consequently, our single-domain spectral 
method cannot be exponentially convergent. Neverthe- 
less, the scheme is rapidly converging towards highly ac- 
curate numerical solutions. 



TO 

ip = 1 + — + u, 
2r 



(28) 



(29) 



(30) 



and we therefore choose to compactify the spatial domain 
by introducing a new radial coordinate A by 



A = [l+ 2r 



(31) 



which implies Ipfjl. 

We separately investigate the two situations in which 
either the linear momentum P or the spin S vanishes. 
A single black hole with either small Bowen-York spin or 
linear momentum has also been considered in I^l40ll41l. 



A. Single puncture with spin 

Consider first a single puncture with P = and 
S l = S x 8\. In the chosen spherical coordinates the 
Hamiltonian constraint reads 



QS" 2 

Aip+ — £r7/;- 7 sin 2 tf = 0. 
I6r b 



(32) 



For the auxiliary function u we obtain a non-linear 
Poisson-like equation, 



UAA 



2u A 1 
~ + ,4 2 (1-A) 2 
36w 2 A(l - A) 2 



Utftf + u$ cot "d + 



sin 2 ■& 



(1 + Au)"' 



sin -d 



(33) 



with w = S x /m 2 . The solution u is uniquely determined 
by regularity and periodicity conditions at i? = 0, d = tt 
and p — 0, <p — 2-7T, respectively. For A = only a 
regularity condition needs to be imposed, while for A = 1 
we set u = 0. Thus, the single-domain spectral method 
described in Sec. [H] is applicable with 

B = 2<&/*-l, (34) 

provided that a global regular solution exists. 
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In order to study the behavior of u globally, and in 
particular close to the puncture, consider the following 
Taylor series which converges for sufficiently small w: 



constraint becomes 



(35) 



All Uj can explicitly be given in closed analytic form. 
In particular, for u\ we obtain (P2 denotes the second 
Legendre polynomial): 



"1 

"1,0 
Ml. 2 



"1,0 + M1,2-P2(C0S#), 

-(-2A 5 + 6A i -5A 3 + 1) 
5 

4 



(1 -A) 3 A 2 



(36) 
(37) 

(38) 



Note that u\ is regular at A — in the spherical coor- 
dinates (A, ip). The same holds for all Uj and in fact 
for u, see |3g. Hence, the C 2 -differentiability of u at the 
puncture has been translated into a C°°-diffcrcntiability 
with respect to spherical coordinates. We can still rec- 
ognize the original behavior of the function which is ex- 
hibited by the fact that u\ possesses a term ~ r 3 which 
is C 2 -differcntiable in Cartesian coordinates. However, 
in the chosen spherical coordinates (A, y), u becomes 
globally C°°. 

Consequently, the application of our single-domain 
spectral method exhibits exponential convergence, which 
can be seen in Fig. ^ for a representative example. 
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FIG. 1: For a single puncture with vanishing linear momen- 
tum parameter the spin S* = m 2 w 8\ with w = 0.2 has been 
chosen. The plot shows the relative global accuracy of the 
spectral method for expansion order ha = nj — n compared 
to a reference solution with n = 50, see (1191 . For this axisym- 
metric example we have used n v = 4. 




qp2 

A^+^rV'" 7 (l + 2cos 2 ^) = 0. (39) 
lor 4 

Similar to the treatment in the previous section we ob- 
tain in spherical coordinates the non-linear Poisson-likc 
equation 



UAA 



2ua 
A 



1 



A 2 {I -A) 2 
q,rM 3 

(l + 2cos 2 tf) 



Utfd + u,g cot 1? + 



sin 2 $ 



4(1 + ,4m) 7 



(40) 



with v = P x /m. Again, we may study the behavior of 
u by performing a Taylor expansion which converges for 
sufficiently small v, 



(41) 



All Uj can explicitly be given in closed analytic form. In 
particular, for u\ we obtain 



Ui 

"1,0 
Ul,2 



"1,0 + "1,2-P2(C0S1?), 

(1-A) 2 



[84(1 - A) log(l - A) 



20A 3 

f 84A - 42A 2 - UA 3 
-7A A -4A 5 -2 A 6 }. 



(42) 
(43) 



(44) 



We recover that u is analytic at A = while it is C 4 - 
differentiable in Cartesian coordinates, which is implied 
by a term ~ r 5 . 

However, the solution u = u(A, <&, ip) also possesses 
logarithmic terms with a branch point at A = 1 (r — » 00). 
For a single puncture, such logarithmic terms are known 
to occur for non- vanishing linear momentum, e.g. [38ll40| . 
In particular, the leading term 



21(1 -A) 3 
5A 3 



log(l - A) 



(45) 



gives rise to a mere C 2 -differentiability of u at A = 1. 
This fact again is reflected by the spectral method, which 
now converges only algebraically to fourth order as ex- 
pected, see Fig. |21 for a representative example. 



TWO-PUNCTURE INITIAL DATA 



B. Single puncture with linear momentum 

Consider now a single puncture with linear momen- 
tum P l = P x 5\ and vanishing spin. The Hamiltonian 



Consider two punctures that are placed symmetrically 
on the x-axis at x = ± b so that D = 2b is the distance 
between the two punctures. We denote the bare mass of 
the punctures by m±, the linear momenta by P± and 
the spin parameters by S±; the subscripts refer to the 
corresponding locations at x = ±6. 
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FIG. 2: For a single puncture with vanishing spin parame- 
ter the linear momentum P r = mv5\ with v = 0.2 has been 
chosen. The plot shows the relative global accuracy of the 
spectral method for expansion order ua = wb = n compared 
to a reference solution with n = 70, see (1191 . For this axisym- 
metric example we have used n v = 4. 



In the following we introduce appropriate coordinates 
in which the auxiliary function u becomes regular at the 
location of the punctures. 

The decomposition l|27|l reads 



2r, 



m_ 
2^ 



u, 



(46) 



with the distances from the punctures given by 

r± = y/[x T b) 2 + y 2 + z 2 . (47) 

As we have seen for the single puncture initial data 
problem, the auxiliary function u discussed there is reg- 
ular at the location of the puncture in spherical coor- 
dinates about this point. We therefore expect a similar 
regular behavior if we were to introduce coordinates that 
become spherical at both punctures. However, regularity 
of u at the punctures can also be achieved if we use spe- 
cific coordinates in which the distances r± are analytic 
functions there (see |23). This is a weaker condition 
because it does not necessarily require one of our coordi- 
nates to behave as r± close to the punctures. 

A coordinate transformation that describes this situa- 
tion at the origin in two dimensions is given by 



where 



C 2 



iy and C = X + iY 



(48) 



(49) 



are complex combinations of Cartesian coordinates (x, y) 
and new coordinates (A, Y) . Clearly, the distance be- 
comes regular with respect to X and Y, 



Note that transformation (|48|l maps a right angle at the 
origin to a straight line through the origin. 

For the two-puncture initial data problem, we apply 
this idea by introducing a specific mapping 



(A,B,ip) i-> (x,y,z), 



(51) 



which is composed of several transformations (see Fig. 

13, 



{A, B,<p) 



(x,y,z). 



(X,R,tp) i-> (x,p,ip) 



(52) 



These transformations are chosen to realize the two dif- 
ferent aspects of the desired entire transformation, (i) 
regularity of r± at both punctures, and (ii) mapping of 
a compact rectangular domain in R 3 to the entire space 
of (x, y, z)-coordinates. 

We first introduce cylindrical coordinates (x, p, cp) such 
that 



y = p cos tp, z — p sin (p, 
and combine x and p to form c, 
c = x + ip. 

Now consider the transformation 



pe [0,2tt), 



C~ 



where C = X + \R. 



(53) 



(54) 



(55) 



It maps the region of the upper half plane with coordi- 
nates {X, R) which is exterior to the unit circle onto the 
upper half plane of our coordinates (x, p), see Fig.^Jc, d). 
The key motivation behind this transformation has been 
to produce locally at each puncture the same effect on an- 
gles as has been done above in the transformation 148|l . 
(|49Jl and which has resulted in the regular expression H50(l 
for the distance from the origin. Similarly we now obtain 
expressions for the distances from either puncture 



r± = \c T b\ = 



2VX 2 + R 2 



R z 



(56) 



\/x 2 + 



y 



'cc 



CC = X 1 



Y 2 



(50) 



which are regular with respect to X and R at the punc- 
tures, i.e. at c = ±6 or C = ±1. 

Next we need to find a transformation which maps a 
compact rectangular region onto the region of (A, R)- 
coordinates. As a first step, the polar transformation 

C = e c , C = C + i^, Ce[0,cx)) ! r 7 e[0 ! 7r] (57) 

yields a strip which is infinitely extended with respect to 
positive £- values, see Fig- EI^ 1 )- Writing c in terms of ( 
gives 

c = 6coshC. (58) 
Thus we recover the transformation 

x = b cosh £ cos r), p = b sinh £ sin 77, (59) 
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FIG. 3: Several coordinate patches for the two puncture initial data problem. Shown are (a) equidistant coordinate lines in the 
system of spectral coordinates (A,B), as well as (b) their images in prolate spheroidal coordinates (£,?;), (c) in the coordinates 
(X,R), and (d) in cylindrical coordinates (x,p). The punctures are indicated by bullets. The (x = 0)-plane, several sections 
of the x-axis and their corresponding images in the other coordinate systems as well as spatial infinity given by A — 1 are 
emphasized by thick lines. 



which maps the well-known prolate spheroidal coordi- 
nates (£, rj) onto cylindrical coordinates. Hence, constant 
£- and fy- values correspond to confocal ellipses and hyper- 
bolas, respectively, in the (x, p)-plane. Their focal points 
are located at the two punctures, i.e. at (£, r/) = (0, 0) and 
(£i?7) = (0i 7r )> see Fig. El The distances r± expressed in 
terms of £ and 77 are 

r ± = fo(cosh £ =f cos 77). (60) 

For a compactification we choose the relations 

7T 

£ = 2 artanh A , 77 = — + 2 arctani?. (61) 

In summary, the transformation from (A, B,ip) to 
(x,y,z) takes the (somewhat symmetric) form 

, A 2 + 1 2B 
x = b ■ 



y = b 



A 2 -I 1 + B 2 ' 
2A 1 - B 2 



I- A 2 l + B 2 
2A 1 - B 2 



I- A 2 l + B 2 



cos ip, 



sin ip. 



(62) 



It is now straightforward to apply our single-domain 
spectral method to solve the Hamiltonian constraint 124|) 
for the two-puncture initial data problem. Again we im- 
pose u — ► as A — * 1, i.e. (x 2 + y 2 + z 2 ) — > 00. As in the 
one-puncture initial data problem, at all the other bound- 
aries we again merely require regularity of the solution 
which replaces a particular boundary condition there. As 
expected, the auxiliary function u is C°° at the two punc- 
tures. 



As mentioned at the end of Section IIIII in general u 
possesses logarithmic terms when expanded at infinity, 
A = 1. In 38] a theorem is proved that does not exclude 
the existence of such logarithmic terms given the fall-off 
condition satisfied by the extrinsic curvature that we con- 
sider here. We have checked in the case of axisymmetry 
analytically that for two equal mass punctures with lin- 
ear momentum logarithmic terms do indeed occur. The 
only exception is when both linear momenta P± vanish, 
in which case the solution is also C°° at A = 1. Otherwise 
we obtain terms ~ (1 — A) 3 log(l — A) if the total momen- 
tum P = P+ + P ^ 0, and terms - (1 - A) 5 log(l - A) 
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if P = 0. In other words, in the center of mass frame 
where the total linear momentum vanishes the leading 
order logarithmic terms cancel, but next to leading order 
logarithmic terms are still present such that the solution 
is only C 4 at A = 1. Although we carried out this anal- 
ysis for axisymmetry, it is to be expected that the same 
result applies to puncture data describing orbiting black 
holes in the center of mass frame. 

A representative convergence rate of our single-domain 
spectral method is displayed in Fig. We show the 
relative accuracy (|19|1 , which involves the maximum over 
a set of points, computed over a 3D set of points as before, 
but also for points only at infinity and at the puncture. 
The error is dominated by errors near the puncture down 
to about 10~ 9 for n < 35. In this regime convergence of 
the maximal error is exponential. However, the error 
at infinity only converges algebraically at roughly sixth- 
order as expected. Around n = 35, the error at infinity 
overtakes the error elsewhere and the overall convergence 
becomes algebraic. 

Therefore we conclude that our numerical method is 
successful since it obtains exponential convergence for 
orbiting punctures down to about 10~ 9 with relatively 
small computational resources. At the punctures our 
coordinate transformation leads to a smooth solution, 
but at infinity there are logarithmic terms which lead 
to algebraic convergence of sixth order. We consider this 
quite satisfactory since higher accuracy is rarely needed 
in numerical relativity. In principle, it should be possible 
to eliminate the leading logarithmic term which should 
bring the calculation close to numerical round-off errors. 
However, it is unclear whether logarithmic terms can be 
avoided completely in this approach, for example by an 
appropriate coordinate transformation. 

In Fig. we compare the result for u obtained by the 
spectral method with u computed by the second order 
finite difference multigrid method on a fixed mesh re- 
finement implemented in BAM previously 0, [^. As 
an example we picked the parameters b = 3M, m + = 
m- = 0.5M, and Pj = -PI = 0.2MS l 2 . The ADM lin- 
ear momentum at infinity vanishes. For the purpose of 
this discussion we have defined M = m++m_. For these 
parameters we can restrict the computational domain to 
one quadrant of a Cartesian box centered at the origin. 
We computed the multigrid data at three overall resolu- 
tions using 7 levels of refinement with approximately the 
same geometrical layout of the boxes. The highest reso- 
lution was obtained for 98 x 98 x 50 points on the finest 
level, while to test convergence we successively doubled 
the grid spacing, resulting in a grid spacing of M/64, 
M/32, and M/16 on the finest level. The face of the 
outermost box is located at about 48M in each case. 

The main result is that the two methods verify each 
other quite accurately on this scale near the punctures 
and also for large x. Near the punctures it is a question 
of resolution whether the known feature of a local inden- 
tation is fully resolved. For the chosen parameters both 
methods reach this level of resolution, but the spectral 
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FIG. 4: Two punctures with vanishing spins. The physical 
parameters are given by m+ = m_ = b, P± — ±0.2 b8\- For 
this plot we took ua = tib = 2n^ = n and compared to a ref- 
erence solution with n — 70. Apart from the global relative 
accuracy (see 11911 ) taken over 6 3 spatial points, the corre- 
sponding maximal deviations at infinity and at the punctures 
are shown. For small n, the error near the punctures is about 
ten times larger than the error at infinity, and the conver- 
gence rate is approximately exponential down to about 10 -9 . 
The error at infinity converges at roughly sixth algebraic or- 
der as expected, and for sufficiently large n this becomes the 
dominant convergence rate. 

method uses significantly less resources. 



VI. TWO PUNCTURES IN THE TEST MASS 
LIMIT 

Apart from the high accuracy that can be achieved by 
spectral methods, they also prove to be very useful for 
the investigation of critical and limiting situations. For 
the binary black hole initial data problem, a situation 
of this kind is encountered when the two gravitational 
sources possess very different masses. It is the aim of 
this section to apply our spectral method for the binary 
puncture initial data problem in this limiting case. 

As a first step we perform the test mass limit analyt- 
ically. The results arising from this study will then be 
compared to those obtained by the spectral method for 
a small mass ratio. 

We consider two non-spinning punctures with bare 
masses m_ and m + located on the y-axis at y = and 
y = D respectively and perform the test mass limit by 
choosing to_ — > with D and m + held fixed. We more- 
over assume that the linear momenta are given by (v- 
and fixed, w_ ^ 0) 

PI = m-v-6i, P| = to_P|, (63) 

which implies that the total linear momentum vanishes 
as m_ — > 0. Thus, in this limit we will have placed 
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(x = 3M, 2 = 0) (x = 3M, 2 = 0) 

FIG. 5: Example for a solution to the Hamiltonian constraint obtained with the spectral method and with a multigrid method 
on nested Cartesian grids. Shown is the regular part u of the conformal factor for two punctures without spin and vanishing 
total linear momentum, which are located on the x-axis at x — ±3Af. Results from the multigrid method are indicated by lines 
with markers. The panels on the left show the various levels of refinement combined into one line for the highest resolution 
(see text). The panels on the right show an enlargement of the region near one of the punctures for three resolutions of the 
multigrid method. In all panels the result for the single-domain spectral method with ua = ub = 40 and n v = 20 is shown as 
a solid line without markers. Note that on this scale the methods agree well both far away and close to the punctures. 



ourselves in a frame in which the total linear momentum 
vanishes. 

In order to understand the behavior of u in the entire 
spatial domain, we have to consider two different ways of 
performing this limit separately: 

1. If we calculate the auxiliary function u at a given 
spatial point at some finite distance from the origin, 
i.e. at fixed coordinates (x,y,z) ^ (0,0,0), we will 
find that u tends to zero in this limit. In particular, 

lim (JL) = *£2L, (64) 
where r — \J x 2 + y 2 + z 2 . The physical meaning of 



the constant A/Zoo (which is obtained through the 
second limiting process, see below) with respect to 
the system's relative binding energy in this limit 
will be discussed in Sec. IVIII 

2. If, on the contrary, we hold the relative coordinates 
(x,y, z) = (x/m-,y/m-, z/mJ) fixed, we main- 
tain finite values for u at these spatial points. In 
particular, the resulting u obeys a constraint equa- 
tion valid for a modified single-puncture initial data 
problem with non-vanishing linear momentum P 
(and no spin). The above constant A/Xoo can be 
read off from these data. 

For both ways of establishing the test mass limit, we 
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rewrite the Hamiltonian constraint as an integral equa- this integral equation becomes 
tion, 

Introducing spherical coordinates 
x = rcosd, 

y = r sin ■& cos <p, (66) 
z = r sin $ sin y?, 

I 



2 2 

m_v_ 
32tt 



2tt 



2r, 



+ 



m_ 
27 



2 dr' 



\x — x ' 



r 



(67) 



where a:' = (x',y',z') with 

x' = r'sini?', 

3SI?'. 



y = r cos c cos <^ , 
z' = r' cos i?' sin ip' , 



and 



r+ = v/.T'2 + (y' - .0)2 + z 
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(68) 
(69) 

and D. They 



The functions g and /i depend on PI, 

remain finite everywhere. 

We now perform the two different limits: 

1. Consider fixed values r > 0. We split the integration 

with respect to r' such that (a) r' e [r/2, oo) and (b) 

r'e[0,r/2]. 

For (a) observe that for r' > r/2 the term 



m + 7Ti— 
2^7 + 27 



remains regular in the limit m_ — > 0, and thus, the con- 
tribution of the corresponding Poisson integral, evaluated 
for r' e [r/2,oo), is of order ©(m 2 ,). 

Performing for the remaining near-zone integral (b) the 
substitution r' = m_s' leads us to 



7rr 



/ V/ 

Jo Jo 



18m_w 2 Z" 271 

' d<p' I 

1 + 2 cos 2 i?' 



7 (2k 



ds' x 



0(s'm- 



[l + 2s'(l + m + /(2D)+u)] 7 



from which it follows that 



lim 



+0 m- 2r 
with the constant A/z^ given by 



(70) 



A^ 



35^2 | J dtpjsmi9(l + 2cos 2 -d)d-d J ds 
oo o 

s 5 [1 + 2s(l + m+l {2D) + u)] -7 . 



Here, the function u is defined by 



lim u(ra_s, ip), 

rri- — >0 



(71) 
(72) 



and turns out to be the auxiliary potential resulting from 
the second limit, which we will discuss now. 

2. Take for the fixed relative distance limit r = m-S 
with s fixed, s > 0, for which we may perform the analo- 
gous steps as in the previous case. We calculate the first 
integral for r' G [£)/2, oo), and again get only a contri- 
bution of order 0(m 2 _). For the near-zone integral we 
obtain 



18i; 2 



d<p' / sintfW 



D/(2m- 



ds' 



X — X' 



, 5 l + 2cos 2 ^ + 0(s / m_) 
[l + 2s'(l + m + /(2L>) + «)] 7 



with vectors 



x = {x 1 y,z) T 
x' = (i\y',S') T 



x/m-, 
x'/rri-, 



where 



x = ssin^, 

y = s cos cos (p, 

z = scostfsmtp, 



x' = s'sini?', 

y' = s' cos -d' cos <p' , 

z' = s'cosi?'siny>'. 
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This leads in the limit m_ — > to an integral equation 
for the function u introduced above. Equivalently, we 
may consider the corresponding differential equation 

9v 2 



Au+ — '(l + 2cos^z?) = 



with 



16s 



■ip = 1 



2D> 



1 

2s 



(73) 



(74) 



and the Laplace operator taken in the spherical coordi- 
nates (s, ip). In particular it follows that 

lim 2su(s, i9, ip) = A/i^. (75) 

S — >OG 

We moreover see that for the function 
rhu with rh = ( 1 



u : 



m+y 
2L>/ 



(76) 



we recover the equation valid for a single-puncture initial 
data problem (without spin), see Sec. IIVBI The (di- 
mensionless) bare mass is just m, and the corresponding 
momentum reads 

TT X = V-lfl 4 . (77) 

The above analytic study shows that we can use our 
spectral methods applied to a single puncture with non- 
vanishing momentum (as described in Sec. HV B|) in order 
to evaluate the test mass limit with algebraic convergence 
of fourth order. These results can be compared with the 
values obtained for a corresponding two puncture initial 
data problem with a small mass ratio, see Table 1. In 
this table one finds the value u_ = u(Q, 0, 0) at the origin 
(i.e. at the 'light' puncture), the expression 



2D 2D , 

u+ = u(0,D,Q) 

m_ 771— 

(i.e. at the 'heavy' puncture), and the limit 

'2ru s 
m,- 



lim 

r — >oc 



(78) 



(79) 



where the latter two tend to as m_ — * 0. We 

have chosen a particular example where the distance D 
and the velocity V- obey the relations valid for the last 
stable circular orbit of a test particle in the gravitational 
field of a Schwarzschild black hole of mass m+ : 



777+ Z 



4V3 
5 + 2V6' 



(80) 



Moreover, we simply set Pi = — v-5\. 

It turns out that for ratios m-/m+ > 1CP 3 the spec- 
tral scheme yields reliable results that approach those of 
the test mass limit. For mass ratios of 10 -3 , four digits of 
accuracy are obtained for the given order of approxima- 
tion from the two-puncture calculation, while six digits 
are obtained with the single-puncture method for the test 
mass limit. 



m-/m+ 




2Du+ /m_ 


lim (2ru/m-) 

r — >oc 


lO" 1 


0.03417 


0.2011 


0.1688 


1(T 2 


0.03406 


0.1635 


0.1601 


10" 3 


0.03406 


0.1596 


0.1592 





0.0340568 


0.159094 


0.159094 



TABLE I: Test mass limit m_ — > for the representative ex- 
ample with values given in 1)80^ with PL — —PL = — m-V-5\. 
For the above non-vanishing mass ratios we used the spec- 
tral method for the binary-puncture initial data problem with 



riA = n B 



2n u 



100. The last line has been calculated with 



the spectral method for the single-puncture initial data prob- 
lem with 7ia = riB ~ 70, n v = 4. 



VII. BINDING ENERGY IN THE TEST MASS 
LIMIT 

In this section we use the results of the previous sec- 
tion to compute the binding energy of two punctures 
without spin in the limit of vanishing mass ratio. The 
aim will be to compare the binding energy in this test 
mass limit with the binding energy of a test particle in 
Schwarzschild spacetime. The deviation of the puncture 
binding energy from the Schwarzschild result will yield a 
quantitative statement about how realistic puncture data 
are in this limit. If punctures were completely realistic 
we should recover the Schwarzschild results. A related 
study of small mass ratios (up to 1/32) has already been 
performed by Pfeiffer (4^ for excision-type initial data 
with Bowen-York extrinsic curvature, and also to a lim- 
ited extent for conformal thin sandwich initial data. 

In order to define a binding energy we need a notion of 
the total mass as well as of the local black hole masses. 
The ADM mass at infinity yields a well-defined global 
mass. For two punctures it is given by 



= 771a 



AMoo, (81) 



where 



1 

'2tt 



Viit dA l = lim 2ru. (82) 



On the other hand, it is impossible to unambiguously 
define local black hole masses in general. In the following 
we choose the ADM mass 



M ADM 



(1 + u±)m z i 



2D 



(83) 



computed in the asymptotically flat region at each punc- 
ture as a measure of the local black hole mass 0. Here 
u + and u- are the values of u at each puncture. As 
shown by Beig [2(| , this definition of local mass has the 
following advantage. For a single slowly moving puncture 
with momentum P^° M = P- , the ADM energy E^ DM 



ADM 



as mea- 



at infinity is related to the ADM mass M. 
sured in the asymptotically flat region near the puncture 
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by 



E 



ADM _ j^ADM 



(Pt DM ? 

2M ADM 



O (P 



,ADM\i 



, (84) 



which is just what one expects if the local mass definition 
is reasonable. If, for example, one uses instead the bare 
mass m_ as the definition of local mass, one finds (e.g. 

Hi) 



rpADM _ m 
tj^ — 171- 



5(P 



ADM\2 



8 m_ 



O 



(iPi DM r), (85) 



which is incompatible with special relativity. 

Next, wc define the binding energy for two punctures 

by 

E b = M ADM - M ADM - M ADM 

m+m_ 

= AMqo — — TO_U_ — — . (86) 

In the test mass limit of m_ — > 0, it follows from Eq. i(64() 
that 



and 



lim AMoo/to_ = A/i c 

m_— »0 



lim u + = 0. 

m._->0 



(87) 



(88) 



Thus Eb goes to zero in this limit. We therefore consider 
Eh/jjL, where 

H = M ADM M ADM / (M ADM + M ADM ) (89) 

is the reduced mass. With the help of Eqs. 1)86(1 . (|64(l . 
l|7H)l and JH3 we find that 



lim — - 



J™o £t \Mf M ' M ADM 
A/ipo (l - g^) ^— 

i i i m_L 

! + «-+ 2D 



(90) 



This binding energy can now be compared with the bind- 
ing energy of a test particle in Schwarzschild spacctimc. 
For circular geodesies in Schwarzschild the binding en- 
ergy, angular momentum and angular velocity observed 
at infinity are given by 



(2r - Mf 



(2r + M)V4r 2 - SMr + M 2 



1. 



L s (2r + M) 2 



Mr 3 



/x 2r 2 V (4r 2 - 8Mr + M 2 ) ' 



and 



Q.c 



' 64Afr 3 
(2r + M) e 



(91) 



(92) 



(93) 



respectively, where M is the mass of the Schwarzschild 
black hole, fi is the mass of the test particle and r is the 
orbital radius in isotropic coordinates. 

In order to compute E^/fx for punctures in the test 
mass limit we have to solve Eq. (|73|l with the appropri- 
ate velocity V- for circular orbits. This raises two ques- 
tions. The first is how to choose the coordinate distance 
D between the two punctures if one wants to compare 
with a test particle in Schwarzschild at isotropic radial 
coordinate r. The answer is that in the limit of m_ — > 
the spacetime is determined by the puncture with bare 
mass 771-1- , so that one simply obtains Schwarzschild in 
isotropic coordinates, which allows us to set 



D = r. 



(94) 



The second question is how one should choose V- for two 
punctures in circular orbit. One could, for example, ob- 
tain V- by requiring equality of Komar and ADM mass, 
which is a necessary condition for the existence of a heli- 
cal Killing vector, as done in An alternative would 
be the effective potential method . Each of these meth- 
ods will give a binding energy and an angular momentum 
for the so-defined circular orbits and in general we do not 
expect the binding energy and angular momentum to ex- 
actly agree with the Schwarzschild results. For simplicity 
and in order to eliminate possible errors in the angular 
momentum we choose 



As A" 



(95) 



so that the angular momentum of the light puncture ex- 
actly equals the angular momentum of a test particle in 
Schwarzschild spacetime. 

Using our spectral method with tia = ns = 70 and 
iij, = 4 we have computed the binding energy for punc- 
tures in the test mass limit. The result is plotted in 
Fig. Inversus the angular velocity given in Eq. (|93|l . Also 
shown are the results for circular orbits in Schwarzschild 
and several other binding energies in the equal mass case 
taken from 01 an d 0]- One can see that the binding 
energy for punctures (solid line on bottom) in the test 
mass limit does not agree with the binding energy of a 
test particle in Schwarzschild (dotted line) , except in the 
Newtonian limit of small Mil. The discrepancy reaches 
about 50% at the innermost stable circular orbit (ISCO) 
of Schwarzschild, which means that the amount of energy 
radiated before reaching the Schwarzschild ISCO is too 
large by 50% and that the location of the ISCO predicted 
by puncture data is wrong. This means that, for the as- 
sumptions made in the definition of the binding energy, 
puncture data are not realistic for extreme mass ratios 
and that one cannot expect to obtain reliable predictions 
about the gravitational waves emitted. 

Let us point out two possible reasons for the discrep- 
ancy. One is that it is not clear whether there are alter- 
natives to our definition of mass, 1)83(1 . that change the 
result. Another issue is that it is known that there is 
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FIG. 6: The solid line on the bottom shows the binding 
energy versus angular velocity for two punctures in circular 
orbit in the test mass limit. For comparison we also show 
the binding energy of a test particle in Schwarzschild (dotted 
line). In addition we show several binding energies for cir- 
cular orbits in the equal mass case. The post-2-Newtonian 
binding energy (broken line) is close to puncture data based 
on an approximate helical Killing vector (pluses) as well as 
to puncture data based on post-2-Newtonian data (squares), 
and also to the Schwarzschild result. 



'artificial' radiation present in puncture data. Such radi- 
ation could contribute at the observed level to the ratio 
of infinitesimal binding energy to infinitesimal mass. 

Interestingly, the curve for puncture data in Fig. [U] in 
the equal mass case (marked by pluses) is much closer to 
both the Schwarzschild (dotted line) case and the post- 2- 
Newtonian (broken line) results, as well as to the results 
of the numerical method based on post-2-Newtonian data 
(marked by squares) discussed in |14|. This might indi- 
cate that artificial radiation affects the binding energy 
per reduced mass for comparable mass puncture data less 
than in the test mass limit. 

Note also that in a method has been described in 
which conformally flat black hole data does indeed lead 
to the correct Schwarzschild result for the binding energy 
in the test mass limit. That method is quite different, for 
example u is approximated by zero and the local masses 
entering the binding energy are defined differently. At 
this point it is not clear how to make contact with our 
approach, but this is clearly an important question for 
future research. 



pected to be at most fourth-order algebraically conver- 
gent. However, one can overcome this problem by intro- 
ducing appropriate coordinates in which the solution is 
smooth at the punctures. In particular, our transforma- 
tion maps the entire K 3 onto a single rectangular domain 
with the punctures at the boundary. 

We have demonstrated rapid convergence of our single- 
domain spectral method and obtained highly accurate 
numerical solutions. Moreover, we have provided a com- 
parison to a numerical implementation with finite dif- 
ferences in Cartesian coordinates and found good agree- 
ment. 

While our coordinate transformation renders puncture 
data smooth at the punctures, in general the fall-off of 
the extrinsic curvature appears to imply the existence 
of logarithmic terms such that the solution is only C 4 
at infinity if the total linear momentum vanishes, and 
only C 2 otherwise. This behavior is a consequence of 
the fall-off of the Bowen-York extrinsic curvature and as 
such unrelated to the puncture construction. It is an 
interesting but to our knowledge mostly open question 
which other approaches to construct initial data for black 
holes share or avoid the problem of logarithmic terms at 
infinity. 

As an application of our spectral method for punc- 
tures, we have considered small mass ratios, and the cor- 
responding results approach the test mass limit which 
was obtained through a semi-analytic limiting procedure. 
Finally, we have computed the binding energy of two 
punctures in the test mass limit and compared it to the 
binding energy of a test particle in Schwarzschild space- 
time and to binding energies in the equal mass case. We 
find that in the test mass limit the binding energy per 
mass deviates from the Schwarzschild result by about 
50% at the Schwarzschild ISCO, while the binding energy 
of two punctures in the equal mass case is close to post- 
Newtonian results, if the ADM mass at each puncture is 
used to define the local black hole masses. This should 
be compared with , where by a different method con- 
formally flat black hole data does lead to the proper test 
mass limit. 

The study of specific coordinate transformations might 
also help in reducing the number of domains that are 
used by methods for binary black hole excision data. 
We have started a corresponding investigation based on 
a coordinate transformation that requires two coordi- 
nate patches, and we intend to apply spectral methods. 
Within the analysis of these data we plan among other 
things to go into the matter of possible logarithmic ex- 
pansion terms of the conformal factor in the context of 
binary black hole excision data. 



VIII. CONCLUSION 

In Cartesian coordinates the regular part of the 
conformal factor of puncture initial data is only C 2 - 
differentiable at the punctures. Therefore, a numeri- 
cal implementation based on a spectral method is ex- 
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